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Abstract. We begin the construction and the analysis of nonosdllatory shock capturing 
methods far the approximation of hyperbolic conservation laws. These schemes 
share many desirable prop er ties with total variation diminishing schemes, 
but TVD schemes have at most first order accuracy, in the sense of 
truncation error, at ex t rem a of the solution. In this paper we construct 
a uniformly second order approximation, which is nonosdllatory in the sense 
that the number at extrema of the discrete solution is not increasing in time. 

This is achieved via a nonosdllatory piecewise linear reconstruction of the 
solution from its cell averages, time evolution through an approximate 
solution of the resulting initial value problem, and averaging of this 
approximate solution over each cell. 


AMS-MOS Classification: Primary 65 MIO, Secondary 65 MQ5 

Key Words: Conservation Laws, Finite Difference Scheme, Nonosdllatory, TVD. 
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1. Introduction. In this paper we consider numerical approximations to weak 
solutions of the scalar initial value problem (IVP) 

(1.1a) u, + f(u) x = u, + a(u) u x = 0 

(1.1b) u(x,0) = uo(x) . 

The initial data uq(x) are assumed to be piecewise-smooth functions that 
are either periodic or of compact support. 

Let v* « v/,( x jt tjf x j 3 J h > *n 3 nr > denote a numerical approximation in 
conservation form 

(1.2a) vf* 1 = vf- Uf)+ v2 ~ fj-vd - (£* ’ A 

Here £* is the numerical solution operator; X = t/A, and fj+w, the 
numerical flux is a function of 2k variables 

(l-2b) /y+ 1/2 = f(vj- k * u ...,vj+ k ) 

which is consistent with (1.1a) in the sense that 

(1.2a) /(«,u,...,u) = /(u) . 

We consider the numerical approximation v A (x,/) in (1.2) to be a 
piecewise-constant function 

(1.3) v h (x,t) = vj, xj.y 2 <x< Xj+yi, nr < t S (n + 1)t . 

Accordingly we define its total variation in x to be 

(1.4) 7V(v") = rn»»(-.0) = j IvjVi - vfl . 


- 3 - 


I 


If the total variation of the numerical solution is uniformly bounded in h for 
0 sfsr 

(1.5) mv*(v)) * C * iT(uo) 

then any refinement sequence h - 0, t = 0(A) has a subsequence A ; - 0 so 
that 

(1.6) v kj — > u 
where u is a weak solution of (1.1). 

If all limit solutions (1.6) of the numerical solution (1.2) satisfy an entropy 
condition that implies uniqueness of the IVP (1.1). then the numerical scheme is 
convergent (see e.g [3] y [12]). 

Recently we have introduced the notion of total variation diminishing (TVD) 

i* 

schemes (see [3]), where the approximate solution operator is required to dimin- 
ish the total variation (1.4) of the numerical solution at each time-step 

(1.7) 7V(v J,+1 ) <s TV(v n ) ; 

these schemes trivially satisfy (1.5) with C = 1. Some early work along these 
lines was done by van Leer in [15]. 

TVD schemes are non-osdllatory in the sense that the number of local 
extrem a in the numerical solution is diminishing in time (as is customary we use 
"diminishing" loosely as short for "non-increasing", throughout this paper). 
Moreover, the value of an isolated local maximum may only decrease in time, 
while that of a local minimum may only increase. 

We were able to construct TVD schemes that in the sense of local truncation 


- 4 - 




error are high-order accurate everywhere except at local extrema where they 
necessarily degenerate into first-order accuracy (see [4], [13], [10], [11], [14]). 
The perpetual damping of local extrem a determines the cumulative global error of 
the "high-order TVD schemes" to be Oft) in the norm, Oft? 2 ) in the L 2 
norm and Oft 2 ) in the L\ norm (see [17]). 

In this paper we introduce a larger class of non-oscillatory schemes, in which 
the solution operator is only required to diminish the number of local extrem a in 
the n umerical solution. Unlike TVD schemes, which are a subset of this class, 
non-oscillatory schemes are not required to damp the values of each local 
extremum at every single time-step, but are allowed to occasionally accentuate a 
local ex tr em u m . 

hi a sequence of papers, of which the present paper is the ilrst, we show how 
to construct non-oscillatory schemes that are uniformly high-order accurate (in 
the sense of global error for smooth solutions of (1.1)). In this first paper we 
describe a second-order accurate scheme of this type. 

The fact that the number of local ex tr em a in the numerical solution may only 
diminish in time is sufficient by itself to guarantee that the application of the 
scheme to monotone data results in a monotone function. Thus non-oscillatory 
schemes, like TVD schemes, are monotcnidty preserving. In particular, when 
applied to a step-function, they do not generate spurious oscillations. 

We note that since the number of local extrema in the solution of non- 
oscillatory schemes is bounded by that of the initial data, uniform boundedness of 
its total variation (1.5) follow*, immediately if the maximum norm of the solution 
is shown to be uniformly bounded. 
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2. Design Principle and Overview 

In this section we describe how to construct a non-osdllatory scheme that is 
uniformly second-order accurate. 

Integrating the partial differential equation (1.1a) over the computational cell 
(*/-i/ 2 > x >+1/2 ) x (t„, r B+1 ) we get 

(2.1a) Of* 1 = 57 - X$ + 1 / 2 (u) “ 

where 

( 2 . 1 b) }j+vM * 7 0 ) * » 

and 

( 2 . 1 c) “7 - \ Jj*” * • 

We observe that although (2.1a) is a relation between the cell-averages uj 
and « 7 +1 , the evaluation of the fluxes f]+\fi(u) in ( 2 . 1 b) requires knowledge 
of the solution itself and not its cell-averages. 

As in Godunov’s scheme and its second-order extension by van Leer [16] and 
Colella and Woodward [2], we derive our scheme as a direct approximation to 

(2.1) . We denote by vj the numerical approximation to the cell-averages uj of 
the exact solution in ( 2 . 1 c), and set vf to be the cell-averages of the initial data. 
Given v" = {vj) we compute v " +1 as follows: 

First we reconstruct u(x,/„) out of its approximate cell-averages {vf\ to the 
appropriate accuracy and denote the result by I(x; v"). Next we solve the IVP 

( 2 . 2 ) v, + f(v) x = 0 , v(x, 0 ) = L(x; v") 
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and denote its solution by v(x,r). Finally we obtain vj l+1 by taking cell- 
averages of v(x,t) 

(2.3) v; +I = i / , '* w v(x,T) dx . 

The averaging operator in (2.3) is non-oadllatory, therefore the number of 
local e xtrem a in v 7,+1 (interpreted as a mesh-function or the piecewise-constant 
function (1.3)) does not exceed that of v(x,t). Assuming v(x,t) to oe the exact 
solution of (2.2) implies (since the exact solution operator is TVD) that the 
number of local extrema in v(x,t) is less than or equal to that of 
v(x,0) = L(x ; v 71 ). Therefore if the number of local extrema in L(x: v n ) does 
not exceed that of v 71 , then the resulting scheme is non-oscillatory. 

We conclude that the design of non-oscillatory high order accurate schemes 
essentially boils down to a problem on the, level of approximation of functions: 
Given cell-averages uj of a piecewise-smooth function u(x), reconstruct u(x) 
to a desired accuracy. Prior to studying this problem we tackle another related 
question in approximation of functions, that of constructing a non-oscillatory 
high-order accurate interpolation of piecewise-smooth functions. 

In section 3 we construct a non-oscillatory piecewise-parabolic function 
Q{x\ u) that interpolates a piecewise-smooth function u(x) at the mesh points 

(2.4a) Q{xf, u) = u(xj) 

and satisfies, wherever u(x) is smooth, 

(2.4b) Q(x\ u) = u(x) + 0(h 3 ') 

Q( x ±0;u) = ± «(x) + Of* 2 ) . 


(2.4c) 





In section 4 we make use of this non-osdllatory piecewise- p arabolic interpo- 
lant io design a non-osdllatory reconstruction of a piecewise-smooth function 
from its cell-averages. As in [16], [2], [5], ind [9] we take L(x\ u) to be the fol 
lowing piecewise-linear function 

(2.5a) L(x; u) = uj + Sj(x - xj)Jh for \x - xj\ < h/2 . 


Unlike the above references that present "second-order accurate" TVD 
schemes, we compute the slopes SJh from Q(x; u ) by 


(2.5b) 


S/h 



Q( x j - 0 ; 5 ), 


^ Q(xj + 0; 5) 



Here m(x,y) is the min mod funoion 
(2.6) m(x 


'■*) = |o 


s • min(|x|,|y|) if sgn(x) = sitn(y) = s 

. otherwise 


We show in section 4 that L(x; u) is a proper reconstruction of u(x) in the 
sense that whenever u(x) is smooth 

(2.7a) L(x; u) = u(x) + 0(h?) 

and 


(2.7b) 


£(*; S) = S(x) + 0(A 3 ) . 


Here u(x) - h~ l u(x + y) dy and L(x; u) - 

* tlf z 

= A” 1 L(x + y; u) dy; like fi(x; u), the latter is also a non-osdllatory 
piecewise-parabolic interpolant of u(x), 

(2.7c) 5) = u(xj) . 


- 8 - 




»/ 


We remark that the "second-order accurate" TVD schemes described in the 
above mentioned reference?* use a slope Sj/h in (2.5a) that approximates 
(i dJdx ) u(xj) to 0(h ) , and their loss of second-order accuracy at local e xt rem a 
points is due to lack of smoothness of the coefficient in the 0(h) term at these 
points 1 . This problem is circumvented in the present scheme by taking Sj/h to 
be (2.5b) which is an 0(h*) approximation to (d/dx) u(xj). Unfortunately there 
is a price to pay for this extra accuracy, namely the loss of the TVD property. As 
in TVD schemes 

(2.8) 7V(v" +1 ) * TV(L(‘\ O) , 

however here 

TV(L('\ v")) a TV(v n ) 

and indeed the scheme may occasionally increase the variation of the numerical 
solution. Although we prove that the scheme is non-osdllatory we have not been 
able as yet to complete a proof of uniform boundedness of the total variation of 
the numerical solution; this is due to lack of techniques to verify uniform bound- 
edness of the maximum norm of the numerical solution. 

In section 5 we study the proposed scheme in the constant coefficient case. 

We verify that it is uniformly second-order accurate, examine its behavior at local 
extrema points and get estimates for the possible increase in total variation per 
time-step. 

In this paper where we consider numerical schemes of the form (1.2) that ar- 

1. We repeat that the results of [8] and [11] imply that TVD schemes, no matter how 
they are constructed, must have this loss of accuracy at local e xtrema 
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derived from approximating the relation (2.1), it is only natural to consider trun- 
cation error in the sense of cell-averages, i.e. we say that the scheme (1.2) is 
second-order accurate if 

(2.9) S"* 1 - E» -5" + 0(A 3 ) 

where u" is the cell-average (2.1c) of the exact solution. Since 

(2.10) a(x) - u(x) + 0(h~) 

whenever u(x) is smooth, (2.9) holds also for pointwise values of the solution. 
However, in the cont e x t of 3rd and higher order accurate schemes, this difference 
in definitions of truncation error will be not only conceptual but of practical 
importance as well. 

Up to this point we have assumed that v(x,t) in (2.3) is the exact solution 
to (2.2). The resulting scheme 

(2.11a) 1 = vj - X[/ y+y2 ( v ) “ fj-\niy)\ » 

where f)+\n(y) is (2.1b) applied to v(x,f), 

(2.11b) it+\nk>) - 7 j£ /«*.<)) it . 

is certainly second-order accurate in the sense of (2.9). Starting with the exact 
cell-averages vj - iij in (2.11) we get from (2.7a) that 

(2.12a) v(x,/) = u(x,f + f n ) + 0(h*) for OsfSr 

and consequently 

(2.12b) //. i/j(v) = + OC* 2 ) . 

which implies (2.9) due to the sufficient smoothness of the coefficient in the 
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Oih 1 ) term in (2.12b). 

} a section 6 we replace the exact solution v(x,r) in (2.3) by an approximate 
one, which we denote by v„(x,/) . This approximate solution is conservative. 
TVD, and second-order accurate in the sense of (2.12a). Thus replacing v(x,/) 
in (2.3) by this approximate solution results in a conservative scheme that is non- 
osdllatory and uniformly second order accurate. 

We remark that an alternative approach to the above is to approximate 
/;+ i/ 2 ( v ) in (2.11b) by using a midpoint rule (or trapezoidal rule) for the 
integral and by replacing v(x,r) with a non-oscillating second-order accurate 
approximate one v„(x,r) (see [16] and [2]). The resulting scheme 

(2.13a) vj * 1 = v? - X(/, +1/2 - fj.vd 

(2.13b) Si*\n. - /(' , »(f/+n 2 . f/2)) 

is certainly second-order accurate, and it is -non -oscillatory in the constant coeffi- 
cient case. Since we have not used the cell-averaging (2.3) to derive this scheme, 
we cannot ascertain in general that the resulting scheme is non-oscillatory. 
Nevertheless, our numerical experiments as well as many other experiments in 
the context of TVD schemes (see e.g. [1], [2]) demonstrate that the numerical 
results are non-oscillatory in many (if not all) applications. 

In section 7 we present some numerical e xp e r i ments that compare the present 
scheme with a typical "second-order accurate" TVD scheme. 

3. Nonoadllatory interpolation. 

The oscillatory nature of second order accurate Lax-Wendroff type schemes 
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results from a Gibbs phenomenon associated with high-order interpolation across 
discontinuities. In this section, as a preparatory step towards reigning a nonos- 
dllatory approximation to (1.1.), we construct a non -oscillatory piecewise- 
parabolic interpoiant Q{x\ u) to a piecewise- smooth function u(x) such that 

(3.1a) Q(x,; u) = m(x,) 

(3.1b) fi(x;u) ■ ^ i+v 2 (x;ii), x ( s:xsx <+1 , 

where ft+ 1/2 is a quadratic polynomial, and 

(3.1c) fi(x; u) - u(x) - O^ 3 ), -j- Q(x ± 0; u) - u(x) - Oih 2 ) 
whvTever u(x) is smooth. 

g(x; 11 ) is non-osdllatory in the sense that the number of its local e xtr em a 
does not exceed that of u(x). 

Since 

“) = a,, qi+vfci+V, «) ' %+i 

it can be written in the form 

(3.2a)q / + 1/2 (x;u) » u, + d (+1/2 u(x - x t )/k + J0/+1/2 a (x - x,)(x - Xt+J/h 2 

where 

(3.2b) 4+ i/2 “ = “f+i - a, 

and D { +y 2 u is yet to be determined. 

0.+ 1/2 “ = ^ “)» x,Sx-Sx l+i . 
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We consider as candidates for q im the two quadratic polynomials 
and qt+i, interpolating u(x) at x h x i+1 ) and (x h x ;+1> x,+ 2 ), respec- 

tively, and choose q\+\n to be the one that is least oscillatory in [x t> x (+1 ]. 
Both qj , j = i and j ~ i + 1, can be written is (3.2a) with M = D j u 

where 


(3.2c) = dj+yjU - dj-yp = Uy +1 - 2Uj + Uj-i . 

Since the least oscillatory of q t and can be characterized as the one that 
deviates the least from the line connecting (x^u,) with (x,+ 1> Uf. f J we choose 
D, + 1/2 u in (3.2a) to be 

(3.2d) Di+yz u = m(D,u, D i+1 u) , 

where m(x,y) is the min mod function 


(3.3) 


m(x,y) 


j • min(|x|,|yD.if sgn(x) = sgn(y) = j 
0 ' otherwise. 


If u(x) is smooth in [xy_ lt xj * 1 ], then qj as a quadratic interpolant of u 
satisfies 


(3-4 )q)(x) - u(x ) = 0(A 3 ), -£■ */x) - -£■ u(x) = O^ 2 ), x y _! s: x s: xj + 1 . 

If Z),u • D 4.JU i 0 then qi+\n is either q t or qi+\. Otherwise we set 
D 1 + 1/2 u = 0, but then smoothness of u implies that Dju = 0(h 3 ) and conse- 
quently qi+ 1/2 - qj = 0(h 3 ) for j = i, i + 1. Thus (3.1c) follows from (3.4). 

We turn now to prove that Q(x\ u) is a nonosdllatory interpolant of u, 
i.e. that the number of its local extrema does not exceed that of u. We do so by 
showing a oqe-to-one correspondence between local e xtrem a of Q to those of the 
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mesh function {uy}, the number of which certainly does not exceed that of u(x). 

Q may have a local extremum in either the interior of some interval 
(x„x; + l) or at a mesh point x,. The first case, which will be refered to as 
interior-extremum, occurs when there is a point x*, x, < x* < x i+1 , such that 

Q(x*; *0 = 0, but qn. v 2 f 0 . 

From (3.2a) if follows that Q has an n tenor-extremum . m (x„ x i+1 ) if and 
only if 

(3.5) \Di+V2 “I > 2|4+l/2 «! • 


qj+i /2 = g,+i/2(jf*)» the value of the interior-extremum is then 


(3.6) 


<?f+l/2 = “i “ 


1 _ 

2 


ltu 


4+1/2 u 
Pi+V2 u 



I 


if Dj+ \j 2 u < 0 it is a local maximum ; if D i+1/ 2 u > 0 it is a local minimum. 
Since D i+U2 u = m(D,u, D /+1 u), (3.5) holds if and only if 

(3.7a) Dju Dt+iii> 0 


(3.7b) | Dj u| > 2|4+ ly2 u(, j = i, i + 1 . 

This implies that q^\a has a locai extremum in (x f , x i+1 ) if and only if both q t 
and 4 ' 4-1 also have a local extremum in (x, , x, + 1 ) and of the same kind. Since 
a parabola has at most one local extrem um, it follows then that q t does not have 
a local extremum in (Xf_i, x,) and 4 + 1 does not have one in (x, +1 , x, +2 ). 
Consequently Q is monotone in both (xj« lf x,) and (x <+1 , x^^, but in an 
opposite sense, i.e. 4-1/2 u ‘ 4+3/2 M < 0; the latter implies that u has a local 
e xtr e m um in [x,-, x, + J and that either u,- or u { +\ is a local extremum of the 


mesh function {i«y} (for obvious reasons the case u, = is counted as a 
sipglg -«TiT «»n i iim ) , The above analysis also shows that In tenor -extrema are iso- 
lated, i.e. if Q has an interior-extremum in (x <t x (+1 ), then it is the oniy local 
extremum of Q in (*,_!, Xj+ 2 )« 

We turn now to examine the case that Q has a local extre m u m at a mesh 
point x,; this will be refered to as a mesh-extremum. The above observation 
that interior extrema are isolated excludes the possibility that Q has an interior- 
extremum in either x,) or (xi,x, + i) and consequently Q is monotone in 

t)yw intervals. This implies that u ' ^i+ 1/2 u < 0 and therefore u; is a 
local extremum of the mesh function {uy}. This concludes the proof that Q(x; u) 
is a non -oscillatory interpolant of u. 

We next express the non-osdllatory nature of Q in terms of total variation. 
If |D/+i/2 “I ^ 2| u j then (2.5) imphes that Q is monotone in (x ; , Xj+J. 
Thus 


(3.8a) \Dj+y 2 m| s; 2\d j+ul “I => 2V[x Jt x,^(Q) - 1 dj+w u| . 

If lDr+i /2 M i > 2|dy+i/2 uj then Q has a local extremum in (x Jt xj+{) and 

= ~ u j\ + i“/+i “ 4j^\n\ • 

Using (2.6) we get 

(3.8b) J Dj+vt <4 j > 2(< ij+yi u| => TV [Xj Jj ^(Q) 


= \4j+U7 M l + i^j+1/2 M l 


dj+l/2 14 
D ]+\n u 



We conclude that 


(3.8c) 0 s 7V(fi) - 2 14-1/2 “ I * 2 |D—l/2 «l f - |) 

; l D m +1 y2 “ 2J 

s 7 2 “i • 

* m€M 

The sum in the RHS of (3.8c) is taken over the set of indices M of intervals 
(**, **+i) in which |D m+1/2 “I > 2|4 ,+i/ 2 «l» where Q has interior- 
extrema. 

Next we show that if u(x) is a piecewise-smooth function of bounded varia- 
tion, then 

(3.9) lim 7V(Q(;u)) = 7V(u) . 

A-0 

We observe that in this case the number of intervals in M is finite and is uni- 
formly bounded by the number of local ext rem a in u(x). Hence (3.9) will follow 
if we prove that D m + 1/2 “ - 0 as A - 0 for all m £ M. To accomplish this we 
show that for h sufficiently small M does not include intervals {x m , x^.^) in 
which u(x) is discontinuous. To see that let us examine the case where u(x) 
has a discontinuity at x € (x„ x 1+1 ) . Q early d t +y 2 u approaches the size of 
the jump in u while J*_. y^u approaches zero as /i - 0. Consequently 

(3.10a) p i u/d i ^ y2 u| = |1 - di^i/ 2 u/d i+ ^ u| - 1 as h - 0 

Hence for h sufficiently small 

(3.10b) 2|d, +w u|> |D f u|i lD, +y2 u| 

which implies i g M. 
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4. Non-orfllatnry Reconstruction. 

Let u(x) be a piecewise- smooth function and denote by u(x) its mean over 
(x — h/2jc + h/2), i.e. 

(4.1) “to - j J “to dy . 

We denote itj = u(xj) and refer to these values as cell-averages of u(x). Given 
{£)}, the task in hand is to reconstruct u(x) to Oft 2 ) in a non -oscillatory way; 
denote the approximately reconstructed function by L(x; u). To achieve Oft i 2 ) 
accuracy it is sufficient to consider L(x; u) to be a piecewise-linear function. To 
make L(x; u) a non-osdllatory approximation we use the non-osdllatory piece- 
wise parabolic interpolation Q(x\ u) to compute its slopes as follows: 

(4.2a) L(x; u) = fij + S/x - xj)/h for \x - xj\< j h 

I* 

(4.2b) Sj = h • Q(xj - 0; 5), Q(x t + 0; S) . 

Here m is the min mod function (3.3); m u and D i+1 ^ u are (3.2b) and 
(3. 2d), respectively. 

We note that L(x; u) may be discontinuous at {xj+ud and that 
(4.3a) L(xj; u) = uj . 

To see that wherever u(x) is smooth 
(4.3b) L(x,u) - u(x) = Oft *) 

we observe that in this case 

(4.4a) u(x) = u(x) + Oft 2 ) 
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and therefore it follows from (3.1c) that 

(4.4b) | Sj - ^ i(xj) + 0(h 2) = ■£ u(xj) + Ofh 2 ) . 

Consequently the RHS of (4.2a) can be expanded as 
(4.4c) L(x; u) = Uj + (x — xj) u(xj) + O^h 2 ) 

= u(x) + (ty 2 ) for |x - x ; | < y A 
and thus (4.3b) follows. 

Denote by L(x; u) the mean value of L(x ; 5) in (x - hJ 2, x + hJ2), i.e. 
(4.5a) L(x; 5) = j J £*£ L(y; u) dy . 

Using (4.2a) to evaluate the integral in (4.5a) we find 
(4.5b) L(x; u) = uj + d J+y2 “ 

• (x - xj)/h + (l/2)(S)+i - 5 ; )(x - xj)(x - xj+d/h 2 , 

for xj*x* xy +1 

(4.5c) L(xy; u) = uj . 

Hence L(x; u), like Q(x; u), is a piecewise- parabolic interpolant of u(x). 
Comparing (4.5b) with (3.2) we find that for x^sxs x^ +1 

(4.6a) L(x; u) - G(x; u) = yCSf+i - 5/ - Dj+vi “)(* “ x j)( x “ • 

From (4.4b) we see that 5, = /i u(x,) + 0(h 3 ) (Note that this is true 


also at local extrema points) and therefore 

S J + 1 " Sj * kl ^2 u (*/+l/2) + 0(A 3 ) . 
On the other hand (3.2) shows that 

-> d 2 

D J+ 1/2 u = h ~^I u ( x J+lt) + 0(A 3 ) . 

Therefore 


^ ' * s !/+l ~ Sj- Dj+wu* o(h 3 ) 

which shows that RHS of (4.6.) is 0(A*). Since (2.1c) show, that 

GO; 5) - «0) = 0(a 3 ) 
we conclude from (4.6a)-(4.6b) that 

( 4 6c ) L(x; u) - u(x) = 0(A 3 ) . 

We turn now to prove that L(x; u) is a non-osdllatory approximation to 

“ (X); m certaial y ™P to «•»* ^(x; 5) is a non-«dlla,ory approximation to 
»(x). We shaU do so by showing that 7V ( , j , j _ J ( i (. ; 5 )), the total-variation of 
L (x, u) in [xj t xy+J, which has the value 


(4.7a) TV [xjtX ^(L( ;u» * 2 (|S)| + ty+jD + \d J+V2 u - j(5y + J >+1 )| , 

can also be expressed as 


(4.7b) 


5)) = max(|d, +w 5|, j|D /+ 1/2 i)) . 


Then it follows immediately form (4.7b), (4.3a) and (3.8) that L is monotone in 
[x y .x /+I ] if and only if Q is; consequently L is a non-osdllatory approximation 
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to u(x) in exactly the same sense as Q is to the interpolated function (see sec- 
tion 3). 

Next let us denote 

(4.8a) Sf = * • Q(xj ± 0; i) , 

i.e. 

(4.8b) Sf = dj- v2 “ + y D j-v 2 Sf = d/+i/2 “ ~ y 0/-*>i/2 u . 

and observe that (4.2b) implies that 

(4.8c) i(|*y| + |S, +1 |) = |[lm(5y,j/)| + W^ 1 .S’/* 1 )|] 

^ yd^/l + 15 + 1 1 ) = *2 - y ^;+i/2 “1 + |d/-t-i/2 M + y ty+1/2 “ 

= maz^iA S|. ylty+i /2 ulj . 

We note that if **| ^ 1^2|°j+i/2 “1 t ^ cn 

sg 0 j^-t-i/2 *• ± *2 0 y+i/2 “ ' *gn(d/+i/2 u) ^ 0 
which in turn implies 

sgn(5j) • sgn(dji +1/ 2 u) 2 * 0, sgn(5 ;+ i) • sgn(d) +1/2 u) * 0 

It follows then from (4.8c) that the RHS of (4.7a) is |dy + y 2 £]. This shows that 

(4.9a) | d J+1/2 “1 * y Pj+V2 “I =» TV^^j(L(-; 5)) = u| . 
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To complete the verification of (4.7b) we still have to show that 


(4.9b) ldj +1/1 S| < i 1 |D,+ W 51 = TV^jCtOi 5)) = y |D.., 5] . 

First we observe that 

(4.10) S+ - Si~ = (di+1/2. « - J £( + 1/2 “) “ W — 1/2 “ + Y £(-1/2 “) • 

= £(“ “ y(^/-t-i/2 “ + £(-1/2 5) 

Since (3. 2d) implies 

(“•ID ^51s±(|D,_ VJ 51 + |D, + ,*5D 

we conclude from (4.10) that 

(4.12) ($,+ - 5f) ' sgn(D, u) * 0 . 

We turn now to prove (4.9b). First let us consider the case that Q(x; u) 
has a local maximum in (xy, x^ +1 ), i.e. Dj u < 0, Dj +1 u < 0, and 

\dj+v2 “1 < y l£y+i/2 “1* 

It follows from (4. 12) that 

(4.13a) Sf a Sj" = d)+y 2 u - y £/+i /2 « > 0 

(4.13b) 0 > dy-m/2 u + y Oy+y2 “ = Sj ~+ 1 ^ <5/^1 • 

The relations (4.13) and the definitions (4.8a), (4.2b) imply that 
(4.14a) Sj = Sf - dj+v 2 “ ~ y £/+i/2 “ 
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(4.14b) 


5/ + 1 - SJ+ 1 = d)+ y2 U + y Dj +lr2 U . 

The same analysis shows that (4.14) holds also for the case that Q(x; u) has a 
local minimum in ( x Jt (4.9b) follows immediately from (4.14) and (4.7a). 

We note that since L(x\ u) is continuous at xj 

(4.15) TV(L(-,u)) = 2 7V M/ _j(L(; !)) = 2 max |4, +1/2 5|, | |D /+1/1 S| 

' ' > 

= 2 I dj+v2 S| + ^2^ 1 D m + yi u] -- |4n + l/2 “l| • 

Here M is the set of indices of intervals (x^,^^) in the interior of which L 
(and also Q) has a local extremum. The number of these intervals is finite and 
is bounded by the number of local e xtrem a of u(x). Comparing (4.9) with (3.8) 
we note that 

(4-16) 7V(L(;S))a7y(C(;ii)). 

5. The constant coefficient case. 

In this section we study the constant coefficient case 
(5.1) u { + oUj — 0, a = const . 

The exact solution of the IVP (2.2) is 

(5-2) v(x,f) = L(x - at ; v") . 

Hence our scheme (2.3) is 

(5.3a) vj> + 1 = L ( x ~ ar; v ") * = L(*j ~ at; v") 
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where L is (4.5a). We have shown in section 3 that the number of local 
extrema in L(x\ v n ) does not exceed that of v n . Since v*"’ 1 in (5.3a) is a cell* 

average of L, it follows that the number of local ertrema in v ** 1 does not 

exceed that of v", and consequently the scheme (5.3a) is non -oscillatory. 

Using (4.5b) in (5.3a) we get the following expression for the scheme 

(5.3b) v/ +1 » (£» • v”)j 

* 

VI - H)(J7 - s 1- 1) * « > 0 

” vf - iu<y+iy2 v" + 1/2 |i(l + |i)(£7+i ~ SJ) if a < 0 ’ 

E h denotes the operator form of the finite difference scheme; ji. = \a, the 

CFL-number, is assumed to satisfy 

(5.3c) W*l. 

We turn now to prove that (5.3) is second order accurate in the sense of 
(2.9), i.e. if u”(x) denotes the mean value (4.1) of u(x,/„) then 

(5.4a) u?* 1 - (£* • 5")/ - 0(h 3 ) . 

To show that we observe that in the constant coefficient case (5.1) 
uf* 1 = iT(xj - ar), and by (5.3a) (E h • = L(xj - ar; u rr ). Hence the LHS 

of (5.4a) is nothing but 

(5.4b) u'ixj - ar) - L(xj - ar; u") , 

which is 0(/» 3 ) as a direct consequence of (4.bc). 

Next we study the time-dependence of the total variation and the maximum 
norm of the numerical solution (5.3). In section 2 we have pointed out that 
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(5.5a) 7V(v" +1 ) * 7V(L(-; v")) . 

Using (4.15) and (5.5a) we get the following upper bound on the possible growth 
of the total variation of the numerical solution per time-step 


(5.5b) 7V(v i, * M ) - TV(v n ) 


2 

miM. 


( 


y HVn/jv-l- ku +u2 v-| 


Here M n is the set of indices of intervals (x m , in the interior of 

which L(x; v n ) (and also Q(x; v*)) has a local extremum. The number of these 
intervals is finite and remains uniformly bounded in time by the number of local 
extrema in the initial data. 

Clearly the upper bound (5.5b) is overly pessimistic. It estimates the possi- 
ble increase in variation in the reconstruction step due to replacing the cell- 
averages vj by the piecewise- linear function L(x; v*). It does not take into 
account the possible decrease in variation in the averaging step (2.3), resulting 
from doing just the opposite, i.e. replacing the piecewise- linear function 
L(x - ar; v*) in (5.2) by its cell-averages (5.3a). 

In the following we shall examine the temporal behaviour of the local 
extrema of the numerical solution and its total variation by analysing the explicit 
values of the cell-averages vf* 1 given by (5.3b). To simplify our presentation 
let us assume that a > 0. 

First we note that (4.8b) implies 

(5.6.) I Sj - Sy-,1 s |S>| + |S,*,| s 2 nuxfl^vl, i 1 Dj.vjv-I) . 

Hence 
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(5.6b) K-L 7 V»| a 1 1 D y _ w V”! = h;. w | - 157 - s;-,|/|dr,. w v"| s 2 
Rewriting (5.3) in this case as 

(5.7a) vf* 1 - v; - M- 4-1/2 V" - in *() - nty-i/2 4-1/2 V" 

“ (1 “ »;-l/2j '7 <T ,- - J7 v 7-l 

where 

(5.7b) aj-vi * m. + y m,( 1 - p.fr./-i/2 • 

We see that the CFL condition 0 < u ^ 1 and (5.6b) imply that 
(5.7c) 1 ; 

thus we conclude 

(5.8) 14-1/2 v”! a Y |D,— 1/2 V*| => v7 +1 € (V7-!. 41 ■ 

Relation (5.8) shows that if v* is monotone for J L ^ j ^ J K , i.e. 

vj L * v Ji * i * ' ’ ‘ * v /,» OT v 4 * v Vi ^ 2 v V tbca vn+1 “ mon °- 
tone for Ji + 1 i j i /*, and in the same sense. Relation (5.8) also shows 
that mesh-extrema of v*, i.e. those for which Q has its local extremum at a 
mesh point, are being damped at the /i-th time-step. Namely, 

(5. 9a) I4 ± i/2 v"| a j V>j±vi V|, i?., => max(v7* 1 , 

(5.9b) !4 ±1/J v”| a I \D uin A, v?., a v? s v7 +1 = mm(v7 +1 . v?*,‘) a v,» . 

We turn now to consider interior local extrema of v*, i.e. those for which 
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Q has its local extremum in the interior of some (x ( , x i+1 ). We recall that such 
a.i extremum is characterized by | y,, l < ^2 v"| and that 5" and 

in this case are given by (4.14); therefore SJVi - 5" * D{* ^ v n . From 
(5.3a) and (4.6a) we see that in general 

(5.10a)vJ'+ 1 1 - Q(*<+i “ ar\v") - y |s(l - \t.)(D t + 1/2 v" - 5{Vi * • 

Hence 

(5.10b) \d, +ul v«| < i |D, + w v*| - vtfj 1 = - ar; v”) . 

Relation (5.10b) confirms the second order accuracy of the scheme at local 
extrema. Although it does not necessitate accentuation of the extremal values, as 
vf# in (5.10b) may still be in [vf, vf +1 ], it dors allow vfo 1 to deviate from 
thi* interval by as much as 

(5.10c) \ |D (+W **|(k+vi V"|/|D,*V7 v”| - m\ 

Thus (5.10b) is the essential difference between the present scheme and the 
"second order" TVD schemes. 

A s imilar analysis, which we do not present here, shows that if vf is a 
mesh-extremum then v^ +1 , j m i, i + 1, relates to Q(xj ~ <jt: v n ) in the fol- 
lowing way: 

(5.11a) v/* 1 x Q(xj - ar; v"), j * i, i + 1, if vf is a maximum 

(5.11b) v/ +1 s Q(xj - aT; v’), ; = /,/ + 1, if vf is a minimum . 

From (5.9)-(5.11) we deduce the following relation between the total varia- 
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tion of the numerical solution and that of its piecewise-parabolic interpolant Q: 
(5.12) TV({Q(xj - ar ; v fl )}) * TV{v H + l ) 22 7V(Q(-; v")) . 

The LHS of (5.12) is the total variation of the mesh function 
{Q(xj - ar, v-)}. Relation (5.12) suggests to consider an equivalent definition 
TV of the total variation of the numerical approximation of the form 

7V(v") = 7V(fi(-; v")) 

with the hope that the scheme (5.3) is TVD with respect to this modified defini- 
tion. Unfortunately our numerical experiments have shown that there are 
instances, although rather rare, that TV(v") is increasing with n; the same is 
true for 7V(v n ) = 7V(L(-; v*)). 

As we have mentioned in the introduction, because of the nonosdllatory 
nature of the scheme, uniform total variation boundedness of the numerical solu- 
tion is implied by uniform boundedness of its maximum norm. If we follow a 
particular local maximum of the initial data we see from (5.9)-(5.10a) that it 
actually decreases most of the time, and whenever it does increase (5.10c) and 
(3.10) suggest that it does so by a "small amount" that vanishes with h -> 0. 

Since the initial data is only piecewise-smooth we have not been able as yet to 
rigorize these arguments. 

We remark that our numerical experiments clearly indicate that in a normal 
computational situation the maximum norm of the numerical solution is indeed 
uniformly bounded. We feel that our inability to prove this fact stems only from 
lack of theoretical tools to analyse pointwise regularity of the numerical solution. 
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6. The nonlinear case. In this section we describe an approximate solution 
v B (x, f) of [5] for the IVP (2.2) 


(6.1) v, + /(v) T = 0, v(x.O) = L(x; v") . 

This approximate solution is consistent with the conservation form of the equation 

(6.1) in the sense that the cell-averaging (2.3) results in a scheme in conservation 
form i.e. 

(6.2) vf* 1 = j v„(x,t) dx = v? - “ fj-\rd 

where the numerical flux fj+m is consistent with /(u) in the sense of (1.2c). 
Furth er mo re , the approximate solution operator is TVD 

(6.3) 7V(v n (-; r)) =s TV(v n (-; 0)) = 7V(L(-; v”))for0s^ T 

and thus by the reasoning presented in section 2, the resulting scheme (6.2) is 
non-osdllatory. 

We turn now to outline the derivation of this approximate solution. To sim- 
plify our presentation we ignore entropy considerations and refer the reader to 
future papers for details of appropriate modifications. We assign to the point 
Xj+\n a characteristic speed that corre s ponds to the Ranldne- Hugoniot speed 
aj+x/2 of the two neighbouring cell-averages vf and v? +1 


(6.4a) 


f /W+1 )-M) 

v J+l - vf 


a )+\Jl 


if vf * vf +] 


a(yf) if vf = vf+! 


and denote by a(x) the piecewise- linear interpolant of {aj +1/2 }, i.e. 


(6.4b) a(x) = a} - 1/2 + («/+i/2 ~ *>- 1 / 2 ) ' 0 “ xj-vd/h 


for Xy _ !/2 ^ x s x/ il/2 . 

The approximate solution v„(x,t) is defined by specifying its constancy 
along the approximate characteristics 

(6.5a) x(f) = x 0 + a(xo) • t 

i.e. 

(6.5b) v ll (x(f), t) - v^xo.O) = L(x 0 ; v 1 *) . 

Using (6.5a) and (6.4b) to ex p res s x 0 in terms of x and t 
6.5C) Xq(x,0 = X]—yi + h (X - Xj-viM/iXj+^t) ~ Xj-wi*)] 


for Xj-^t) < x < x j+l/2 (t) 
we get from (6.5b) that 


(6.6a) v„(x,/) 


* 

= L X J- 


in 


+ h 


* ~ 1 / 2(0 
*/+ 1/2(0 ” 3 " 1/2(0 * 



for Xy.j^O) < X < 1/2(0 


where 

(6.6b) *<+i/2(0 = *i+i/2 + r • ai+M . 

Taking cell-averages of the approximate solution (6.6) we get the conservation 
form (6.2) 

(6-7a) v} +: - vj ” ^(/y+i/2 ~ fj-w) 
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with the numerical flux 


(6.7b) 



f(vf) + 1/2 aj +] j 2 (1 - X ajf-i/ 2 ) * Sj 
f( v j+l) ~ 1^2 aj+u2(l + V ■ Sj+i 


v 


if in ^ 0 

if 1/2 s ^ 


where 


(6.7c) i, = 57/(1 + X(S } +1 /2 - ay- 1 / 2 )] • 


Note that (6.7) is identical to (5.30) in the constant coefficient case. 

We turn now to prove that the scheme (6.7) is uniformly second-order accu- 
rate in the sense of (2.9). Starting with the exact cell-averages vf = uj in (6.7) 
this amounts to showing that 

(6.8a) //+ 1/2 = fjrvM + Oik 1 ) 

with a sufficiently smooth coefficient in the OQ 1 2 ) term; here fj+\n is the 
numerical flux (6.7b) computed with the exact cell-averages, and fj+ 1 / 2 OO is 
(2.1b). We shaU do so in two steps: first we shall show that 

(6.8b) = y !/(£(*/+ 1 / 2 ; “")) + f(L( x o( x j+\rt’ T);iT))] + Oih 2 ) , 


where x^xj+^t) is (6.5c), and then we shall verify that 

(6.8c) y[/(L(xj +1/2 ; 5")) + /(L(xq(x, + Wf r); 5"))] = /; +y2 + Oih 1 ) . 

Special attention will be given to the smoothness of the Oih 1 ) coefficients . 

To show (6.8b) we start by using the trapezoidal rule to approximate the 
integral in (2.1b); we get 

(6.9a) fjM j 2 (u) = y [/■(“(*,< +i/ 2 > '*)) + /(“(*y+i/ 2 > *n + t))] + O^ 2 ) . 
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The smoothness of the O^h 2 ) term follows from that of /(a) and u(x,/). Next 
we observe that a(x) in (6 4b) approximates a(u(x,f „ )) to 0(h?), and there- 
fore we can use the approximate characteristic line (6.5c) to trace 
u(x ; + 1 / 2 , t n + t) to u(x Q (x j+l/2 , t), t n ) with 0{h 3 ) accuracy; consequently 

(6.9b) f(u(xj+v2, t n + t)) = f(u(x 0 (xj+v2> t),' Q) + 0(A 3 ) . 

Finally we obtain (6.8b) by approximating u(x,t n ) in (6.9a) and (6.9b) to 
0(h*) by I(x; iT) (see (4.4)). The smoothness of the 0(h*) term in this 
approximation is due to (4.4c): 

SJ = h ■ «,(*, «.) + 0(* 3 ) . 

(We recall that the degeneracy to first order accuracy at local extrema points of 
some "second-order accurate" TVD schemes is due to lack of smoothness there of 
the 0(h *) term in (2.7a)). 

We turn now to verify (6.8c). First let us consider the case oj + w a0: 

O = “7 + \ SJ , 


£(*oC*)>i/ 2> t ); «") = y} + 



\a 


7+U2 


1 + ~ * j - 1/2) 


157 


> 


and expand the LHS of (6.8c) around uj 1 . We get 

(6.10a) /(ip + i o(up(l - X S J+ v,) Sj + jo'(5p[(l - X Sj.^ 2 


+ (Xa^+ia) 2 ]^) 2 + 0(h 3 ) = fi*\n + Y (2X 2 o 2 - 1) • o’ • («,) 2 |/.n /2 + 0(X 3 ) . 


Similarly in the case aj+ 1/2 ^ 0: 
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U*j+ 1/2J «") 


U J + 1 9 


£(*c(*j+i/2> T ); 5") 




Xaj-H/2 

1 + \(aj+y2 *]+\n) J 


SJ + 1 . 


we expand the LHS of (6.8c) around uj+i to get 
(6.10b)/(u7 +1 ) - y a(5/ +1 )(l + Xa) J . 3 ^)iy* 1 

+ J «'(“7+l)[(l + + (Aaj+i/J 2 ] • (5j +l ) 2 + 0(A 3 ) 


- /j+l/J + Y (2AV - 1) • a' • W 2 (,* w + 0(* 3 ) . 

We see from (6.10a) and (6.10b) that independently of the sign of 5) +1/ 2 , the 
0(A 2 ) term in (6.8c) is the same, namely 

« * 

Y ( 2x2flZ “ *) * ^ * (“x) 2 l I/+ 1/2 • 

This completes the proof that the scheme (6.7) is second-order accurate in the 
sense of (2.9) wherever u(x,t ) is smooth, including local extrema and sonic 
(f = 0) points. 

Remarks : (1) The numerical flux (6.7b) can be rewritten as: 

(6.11) fj+ 1/2 = + Avf+l) - F,+wl (v; + i - v7) 


+ max(0, fly+ 1 / 2 ) ’ (1 “ Xay_i/ 2 ) * Jy 


— min(0, 5/+i^) * (1 + Xa^+^t) • 
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(2) Our proof that the scheme (6.7) is non-osdllatory is based on the 
representation of (6.7) as the cell-average (6.2) of the non-osdllatory approxi- 
mate solution v n (jr,x) in (6.6). To ensure that v„(x,t) remains univalued for 
OifSr we have to restrict the time-step t so that for all ; 

(6.12a) ^(t) > xj-v2( t) . 

This implies the condition 

(6.12b) t maXj(aj-v 2 ~ a j+ 1/2) ^ ^ • 

On the other hand, to derive the particular form of the numerical flux (6.7b) we 
have made the assumption 

(6.13a) |xy+ 1 / 2 OO - xj+Lal < h for all ; , 

which implies the condition. 

(6.13b) t max j\aj+ j^l < h . 

The two time-step restrictions (6.12b) and (6.13b) are characteristic to 
Godunov-type schemes. The practice of reconciling the two conditions by 

t maxy|a]i + 1 / 2 | ^ h/2 

is completely unnecessary: The scheme (6.7), as the original Godunov scheme, 
remains non-osdllatory (or TVD in the case analysed in [10]) for the full CFL- 
restriction (6.13b). The reasoning for this observation is as follows: The approxi- 
mate solution (6.6) under the CFL restriction (6.13b) may fail to be univalued in 
they-th cell only if a)- 1/2 > 0 and a}* 1/2 < 0. In this case the numerical fluxes 

A 

fj ± 1/2 as defined by cell- averaging in the neighboring cells j ± 1, remain (6.7b). 
Thus the only thing that needs to be changed in this case is the definition of 


v„(x,t) in the y'-th cell. 

(3) We observe that once v„(x.t) is defined globally in (6.6) there is no 
intrinsic need to average it on the original mesh. We may average it on different 
intervals and still conclude that the resulting approximation is non-osdUatory and 
conservative. Furthermore, the construction of the interpolant Q, the approxi- 
mation L and the approximate characteristic field a(x) needed to define 

v„(x, r), does not depend on the uniformity of the mesh. Therefore the scheme 
(6.7) generalizes immediately to non-uniform moving meshes. Of particular com- 
putational interest are the self-adjusting moving grids of the type described in 
[12], which make it possible to obtain perfectly resolved shocks and contact 
discontinuities. 

(4) We note that since the approximate solution v„(x,t) in (6.6) is conser- 
vative, it is possible to consider an associated random-choice method obtained by 
replacing the cell-averaging in (6.2) by a sampling with respect to a variable that 
is uniformly distributed in the cell, i.e. 

»7 +1 = V.(xj + 8/A,t) 

where Sj 1 is uniformly distributed in [-1/2, 1/2]. Following the reasoning of [7] 
it is clear that the resulting random-choice method is non-osdUatory and that its 
limi ts are weak solutions of (1.1). Although the random-choice approach has 
many attractive computational features, it has been our experience that in many 
application it is possible to accomplish the same computational goals with a self- 
adjusting moving grid. In this case the use of the latter is preferable as it offers 
gain in resolution without a loss in pointwise accuracy that is associated wiJ< sam- 
pling. 
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7. Numerical Illustration. In this section we compare the new uniformly 
second-order non-or»dllatory scheme of this paper (to be refered to as UN02) to 
the typical second-order TVD scheme (to be refered to as TVD2). Both schemes 
can be written in the form (6.7), i.e. 

(7.1a) yT* “ Y? - ^lW^-i/s) , 


(7. lb) 1/2 


fiyj) + 1/2 5 )+ 1/2(1 ~ i/ 2 ) 5 "/tl + X (a^ 4. 1/2 - i/2)l 

tf 1/2 ^ 0 

f(?j+ 1) ~ 1/2 fljf+l^(l + ^+3^)5"+ i/(l + H<*j+y2 ~ <*J+V2)] 

5/+1/2 S 0 


(7.1c) 57 - m(5/^“); 

here ay+ 1/2 is (6.4a) and m(x,y) is the min mod function (3.3). Sf are dif- 
ferent for TVD2 and UN02: 

(7.2) TVD2: Sf =<< /± 1 / 2 v" 


(7.3) UN02: Sf = =F \ D JiUi y" , 

where 1/2 and are defined in (3.2). 

UN 02 and TVD2 are both second-order accurate Godunov-type schemes 
that differ only in the reconstruction step (4.2a): 

(7.4a) L(xyi) - u< f Sj(x - xj)/h for \x - xj\ < h/2 , 

when’; the slopes of the lines are calculated by (7.3) and (7.2), respectively. 
Therefore we start our comparison on the approximation level. 

In Table 1 and Fig. 1 we present approximations to 
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u(x) * siniTx, -lsxsl. Wc divide [-1,1] into V equal intervals and 
define 

(7.4b) j:. 0 sjxN. 

The symbols in Fig. 1 denotes values of uj - simrxy for N * 10 in (7.4b). In 
Fig. la we show the piecewise-parabolic interpolant Q(x jt) (see section 3). In 
Fig. lb we show the piecewise-linear approximation L^x^i) which is (7.4a) 
with (7.1c) and (7.3). In Fig. lc we show the piecewise-linear approximation 
lTVD^xju) which is (7.4a) with (7.1c) and (7.2). We make the following obser- 
vations regarding Fig. 1: (i) Q is a better approximation than I UN02 ; I UN02 is a 
better approximation than L TVD2 . (ii) 7V(L UN02 ) > 7V(u) > 7Y(L TVD2 ). hi 
table 1 we quantify the first observation; we list the L«-e rror and the 12 -error of 
these approximations to simrx for a refinement sequence of N = 10,20,40,80 
in (7.4b). Gearly Q is an 0(h 2 ) approximation, while L UN02 and L TVD2 are 
i 2 ). The error in L UN02 is about a 1/3 of the error in L TVD2 . 

In Table 2 and Fig. 2 we present solutions of UN02 and TVD2 for the con- 
stant coefficient case 

(7.6) u, + u, * 0, u(x,0) = sin irx, -lsxs 1 

with periodic boundary conditions, at t = 2 with r /h - 0.8. Figs 3a and 3b 
show UN02 and TVD2, respectively, with N =* 20 in (7.4b). In table 2 we list 
the I.-error and Lj -error for a refinement sequence with N - 20,40,80,160. 
Gearly UN02 is second-order accurate in both L m and L lt while TVD2 is 
second-order accurate in L\ but only first order accurate in L*. Fig. 3b demon- 
strates that the degeneracy to first order accuracy at local extrema in the TVD 
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scheme adversely affects the accuracy everywhere (Because the scheme is TVD it 
cann ot switch ab. aptly to second-orr* t accuracy as this would introduce osdlla- 
tious; consequently it takes quite a while to recover the second order accuracy). . 

Next we approximate the discontinuous function 


(7.7) 


u(x) 


-xsin(3irx 2 /2) , -1 < x < -1/3 

|sin(2irx)| , |x| < 1/3 

2x - 1 - l/6sin(3irx) , 1/3 < x < 1 


which we extend to have period 2 outside [-1,11- 

In Fig 3 we present approximations to <J>(x), using N = 20. Fig 3a shows 
Q(x;u), Fig 3b shows L UN02 (x;5), and Fig 3c shows ^ TVD 2 (x;m). We again 
observe that Q is better approximation then L UN02 , while L 1 "™ 2 is a better 
approximation then L TVD2 . 

In Fig 4 we present solutions of UN02 and TVD2 for the constant coeffi- 
cient problem (7.6), initial data given by (7.2), and periodic boundary conditions. 
We take t = 2 and j/h = 0.8. Figs 4a and 4b show UN02 and TVD2 respec- 
tively with N = 40. Fig 4b shows the damping effect that the TVD scheme 
imposes due to its degeneracy to first order accuracy at local extrema. 

In Fig 5 we solve the same problems, except we impose boundary conditions. 
At x = - 1 we impose the given function (7.7) evaluated at - 1 - r. No boun- 
dary conditions are imposed at x = 1. We implement this numerically using 
UN02 and TVD2 except at the boundary points. There we are in general, unable 
to construct non-osdllatory piecewise parabolic interpolants Q(x,u), so we con- 
struct the only possible parabolic intr.polant thru xj,x <+1 and the point to either 
the left or right which lies in the region. The analogous procedure is carried out 
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at the reconstruction stage. Figs 5a and 5b again show the results at t - 2 with 
r/h » 0.8. 


The possible introduction of oscillations through the boundary conditions 
does not seem to have degraded the p e rformance of either scheme (in fact the 
opposite is observed). Again the TVD2 scheme shows a damping effect. 

In Table 3 and Fig. 6 we present results for Burgers’ equation 


(7.8) u t + utij - 0, u(x,0) - a + sin ^(x +0), -lSxSl 


with periodic boundary conditions and r Zh (1 + |a|) - 0.5. The solution to (7.7) 
is smooth for t < 1/it; at t » 1/tt it develops shocks. In Table 3a we list the 
La -error and Lj -error of UN02 and TVD2 at t =* 0.15 for a = 0 = 0 in 
(7.7). This table shows the same asymptotic behaviour as Table 2, except that 
because of the large gradients it shows for a smaller h. 

In Figs 6a and 6b we show results of UN02 and TVD2 for (7.8) with a = 2 
and 3 * 1 at / « 0.35 with N =* 20. In this case the solution to (7.8) 
develops a shock moving with speed 2 beginning at time r » 1/ir = 0.318. 

In Table 3b and Figs 6c and 6d we repeat the previous calculations for the 
schemes (2.13): 

( 7 . 9 ») V?* 1 -vj- *(jl+V2 - fj-U2) 


(7.9b) /j+i/j - /(»„(*;+ 1/2, T/2)) - f(L Mij+m, t/2),v")) 

where x^xy+yj* t/ 2) i* (6.5c), i.e. 


(7.9c Yj+vi = 


/ 


v7 + 



1_ 1 ~ X(dy+j/2 + 

2 1 + X(dy +1/2 ~ &j-yi/2 { 

l_ 1 + H&j+in + n/2 

2 1 + \.(dj+i/2 ~ &j+viy2 J \ 


if <*y+i/2 ^ 0 
if dy* 1/2 ^ 0 
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As we have remarked in section 2. v }* 1 in (7.9a) is not a cell-average of 
v „(jr,r), but only an approximation to it. Therefor* it is not necessary to take 
dj+v2 (7.8c; to be (6.4a). We choose so that (7.9c) is continuous at 

=* 0 : 


(7.9d)Vy2 


l/tf,; - 7- c ;.i) - f(y] + \sf) 


- }sr*i) - (*/ + \sj) 


We denote the schemes (7.9) with 5" defined by (7.1c) and either (7.2) or 
(7.3) by FVD2 and FN02. respectively. We note that (7.9) is identical to (7.1) in 
the constant coefficient crje, and consequently FVD2 and FN02 are nonoscilla- 
tory in the constant coefficient case. Figs 6c and 6d show that FN02 and FVD2 
are also non oscillatory in the case (7.8). Furthermore, Table 3b shows that 
FN02 is much more accurate than UN02 (FVD2 is about the same as TVD2). 

In all previous examples we have pr es ented pointwise calculations; namely, 
we have initialized the numerical solution by taking vP to be the value of the ini- 
tial data at x ; , and we have considered vj to be an approximation to u(x ; ,r n ). 
(Surely this is an acceptable practice for second order accurate schemes.) In Table 
3c we repeat the calculation for UN02 in Table 3a, but now in a sense of cell- 
averages and denote it by AN02. Now we initialize UN02 for (7.8) with 
a = 0 = 0 by cdl-averages of the initial data, i.e. 

(7.10a)vj ) - (cos (vxj+vd - cos^xy.^] = • *in(irx y ) , 

and regard vj to re p resent cell-averages of u(x,f„). To obtain a pointwise 
approximation to u(x,r„) we first compute point values ^ of its indefinite 

integral u m (xj+yd - f ^“CV*) <*y by 

Xq 
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Next we obtain a global piecewise-linear approximation v(x,/„) to u(x,t n ) by 
(7.10c) »(*,/„) = ■£ e(*;V") 

where Q is the piecewise- parabolic interpolant of section 3. Finally we get 
(7.10d) Q(XJ-,V") = i (v; +1/ j - '/;. w ) = yj . 

Tims the only difference between AN02 in Table 3c, and UN02 in Table 3a 
is the initialization (7.10a), which itself differs only slightly from the mesh values 
of the initial data (since sin(irA/2)/(irA/2) = 1 - l/6(ir/i/2) 2 + 0(h 4 )). 

We remark that cell-averages do play a significant role when the initial data 

is discontinuous (since they provide information about tne location of the discon- 

, ♦ 

tinuity) and in higher-order Godunov-type schemes; this will be described else- 
where. 

We turn now to present calculations with a formal extension of UN02 and 
TVD2 for systems of conservation laws. We consider a Riemann problem for the 
Euler equations of gas dynamics 



(u L x < 0 

(T.ila) 

“f + Mx = “(*»&) = L r x > o » 

(7.11b) 

u = (p,m,£) r , /(u) = (m,m 2 /p + P, m(E + P)J p) r , 

(7.11c) 

p = (7 - 1 )(£ - \ p) . 


Here p, m, E and P are the density, momentum, total energy and pressure, 


I 


respectively; we take y = 1.4. 

In the following we apply the extension technique of [3] to UN02 and TVD2. 
The idea is to extend UN02 and TVD2 to systems in such a way that will be 
identical to (7.1) in the scalar case, and will decouple into (5.3) for each of the 
characteristic variables in the constant coefficient system case. To accomplish 
that we use Roe’s averaging for (7.11) (see [13]) 

(7.12a) Vj 4 .y 2 * V ( v j> v 7+i) 

for which 

(7.12b) /(yj+i) “ f(yj) = a ( v j+i/i)( v 7+i “ v j)> A(u) = d f /du , 

and define local characteristic variables with respect to the right-eigenvector sys- 
tem of i4(vy +1/2 )- We extend ( 6 . 11 ) to systems as follows: 

(7.13a) vf +l = v/ - X(fc +ly2 - fj-vJ 

(7.13b) /, + 1/2 = | j/(>f) + /(•?,,) - i cf +y2 Rf +l/1 

(7. 13c) — max(0.a , _. 1/2 )(l — \aj^ 1 / 2 ) 5 ? 

+ j^)(l + \aj+yi)Sj^i . 

Here aj+in is the i-th eigenvalue of 4(v ;+1 ^) corresponding to Rj+y-i, and 
dj+ i/2 w denotes the component of dj+y^v = v ™ +1 - vj 1 in the i-th characteris- 
tic field, i.e. 

(7.13d) dj+y2 v = 2 (^*+ 1 / 2 h ')^;> 1/2 • 
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Likewise Sj denotes the component of the vector of slopes in the ifc-th charac- 
teristic field, and is defined as follows: 

(7.13e) Sj = m(S l Lj,S+j)/[ 1 + k(<#+i /2 - <^-\n)\ I 

m(x,y) is the min mod function (3.3). S±j are different for VD2 and N02: 

(7.14) TVD2: S^j = df ±u2 w 

(7.15) UN02: S k ±J = d} +vl w =F ± D} ±v2 w, 

Di±V2 w = m(df+y2 w - df+viw, df+\n w ~ <#- 1 / 2*0 • 

In Figs 7.8 and 7.9 we show numerical solutions of UN02 and TVD2, respec- 
tively, for the Riemann problem (7.1b) with 

U L = (l,0,2.5) r , U R = (0.125,0,0.25) . 

it 

These figures demonstrate that the formal' extension to systems is nonoscillatory 
in this case. Since the solution to the Ri emann problem is just constant states 
seperated by waves we do not get to see here the extra resolution power of UN02, 
except that its numerical solution is somewhat "crisper" than that of TVD2. In 
this calculation we have not employed any artificial compression in the linearly 
degenerate field and therefore the contact discontinuity smears like n 173 , as 
expected. The interested reader is referred to [4], [5] and [10] for a detailed 
description of such compr es sion techniques, as well as for details of entropy 
enforcement mechanisms. 
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Table 1 : Approximations to u(x) = sin(iTx), — 1 as x ss 1 , with periodic boundary conditions. 



Loo-ERROR 

L r ERROR 

N 

Q 

^UN02 

^TVD2 

Q 

^UN02 

L TVD2 

10 

1.545 x 10“ 2 

5.122 x 10“ 2 

1.420 x 10" 1 

1.494 X 10“ 2 

2.467 x 10“ 2 

7.016 x 10" 2 

20 

1.971 x 10" 3 

1.231 x 10 -2 

3.558 x 10" 2 

1.802 X 10“ 3 

5.576 x 10" 3 

1.525 x 10" 2 

40 

2.476 x 10" 4 

3.083 x 10“ 3 

9.163 x 10~ 3 

2.148 x 10“ 4 

1.355 x 10“ 3 

3.902 x 10“ 3 

80 

3.104 X 10“ 3 

7.710 X 10“ 4 

2.308 x 10“ 3 

2.617 x 10“ 3 

3.351 x 10“ 4 

9.787 x 10“ 4 


Table 2: Numerical solutions of u, + u, = 0, u(x,0) = sin nx, -1 s x s 1 at t = 2 
with periodic boundary conditions and r/h = 0.8. 



Loo-ERROR 

L r ERROR 

N 

UN02 

TVD2 

' UN02 

TVD2 

20 

7.097 X 10“ 3 

8.119 x 10“ 2 

8.944 x 10" 3 

6.778 x 10“ 2 

40 

1.607 x 10" 3 

3.477 x 10“ 2 

2.044 x 10“ 3 

2.033 x 10" 2 

80 

3.870 x 10" 4 

1.453 x 10“ 2 

4.926 x 10“ 4 

5.626 x 10“ 3 

160 

9.201 x 10” 5 

5.975 x 10" 3 

1.172 x 10“ 4 

1.528 x 10“ 3 




Table 3a. Numerical solutions of u, + uu x = 0, u(x,0) = simrx, at t - 0.15 and r/h - 0.5 - 
UN02 and TVD2 



I.-ERROR 

L r ERROR 

N 

UN02 • 

TVD2 

UN02 

TVD2 

"20" 





40 

5.712 x 10" 3 

1.054 x 10" 2 


5.051 x 10’ 3 


1.552 x 10“ 3 

4.422 x 10 -3 


1.340 x 10" 3 

160 

3.985 x 10" 4 

1.837 x 10“ 3 

1.965 x 10“ 4 

3.621 x 10“ 4 


Table 3b. Same as Table 3a for FN02 and FVD2. 



L.-ERROR 

Lj-ERROR 

N 

FN02 

FVD2 

PN02 

FVD2 

20 

| 




40 

1.959 x 10" 3 


8.869 X 10"“ 


B 

5.106 x 10“ 4 


2.163 x 10" 4 

1.072 x 10“ 3 

160 

1.251 x 10” 4 



2.946 x lO -4 
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Figure la 

~i-rure 1 . .'.DDr-'xitr.atirns ?f u = sinnx. -1 ^ x * 1, with N - 10. 
a) Piecewise-Parabolic interpolant Q(\,u). 



Piecewise-linear approximation L TVD2 (x*u) 
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FlgureJ^. Numerical solution of u,. + U)[ - 0, u(x,0) Figure lib . Same as Figure 4a using TVD2. 

defined by (7.7) at t 2 with periodic 
boundary conditions N = ho f/h - 0.8, 
uainR UN02. 










Figure 6a. Solution to (7*8) with 01 
CFlJ =.y, at t = 0.35 using UN02. 







Figure 6C, Same as Fig. 6a using FN02. Figure 6P . Same as Fig. 6a using 




£i£- ure Figure 7b . 

Figure 7 . Numerical solution of density in a RLemann problem for Euler's equations, 
(a) UN02 (b) TVD2 
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?*jc-.erical solution of velocity in a Rtemann problem for Ruler's equations. (a) UIJ02. (b) T D2. 












